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ABSTRACT 

Many astrophysical sources of high energy emission, such as black hole magneto- 
spheres, superstrongly magnetized neutron stars (magnetars), and probably relativistic 
jets in Active Galactic Nuclei and Gamma Ray Bursts involve relativistically magnet- 
ically dominated plasma. In such plasma the energy density of magnetic field greatly 
exceeds the thermal and the rest mass energy density of particles. Therefore the mag- 
netic field is the main reservoir of energy and its dissipation may power the bursting 
emission from these sources, in close analogy to Solar flares. One of the principal 
dissipative instabilities that may lead to release of magnetic energy is the tearing in- 
stability. In this paper we study, both analytically and numerically, the development of 
tearing instability in relativistically magnetically-dominated plasma using the frame- 
work of resistive magnetodynamics. We confirm and elucidate the previously obtained 
result on the growth rate of the tearing mode: the shortest growth time is the same 
as in the case of classical non-relativistic MHD, namely r — ^TaTd where is the 
Alfven crossing time and t^^ is the resistive time of a current layer. 
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■ 1 INTRODUCTION 



Dissipation of magnetic energy may power high energy 
■ emission in a variety of relativistic astrophysical phenom- 
' ena, e.g. pulsar wind nebulae llUoroniti 19901 lUsov 19941 
[Lyubarsky fc Kirk 200T| |Kirk fc Skjaraasen 200 3|, jets 
of active galactic nuclei jRomanova & Lovel ace 19921 
[Jaroschck ct al. 2001), gamma- ray bursts il'Drenkhahn 20021 

^renkhahi^&^Sprui^JLC^^^ and magnetars 

^Thompson fc Duncan 1995| [Lyutikov 2003| l. 

Perhaps the most clear-cut case for magnetic dissipa- 
tion is magnetars - neutron stars with superstrong magnetic 
fields , sometimes in excess of quantum magnetic field Bq = 
mlc^ /he = 4 X lO^^G. (Observationally, magnetars are sep- 
arated in two classes - Anomalous X-ray Pulsars (AXPs) 
and the Soft Gamma-ray Repeaters (SGRs) - both showing 
X-ray flares and quiescent X-ray emission.) A number of ev- 
idence suggest that processes that lead to the production of 
magnetars X-ray flares (and possibly of the persistent emis- 
sion) are similar to those operating in the Solar corona. The 
bursting activity of SGRs is strongly intermittent, showing 
a power law dependence of the number of flares on their en- 
ergy, dN/dE ~ E" , with a = 1.66 and log- normal distribu- 
tion of waiting times between the flares l |Gogiis et al. 1999l l, 
both being similar to Solar flares. Giant flares (GFs be- 



black hole physics - MHD ~ methods:analytical - meth- 



low), immense explosions releasing up to 10*^ ergs in a 
fraction of a second, provide a number of independent ev- 
idence in favour of magnetic dissipation. Two recently ob- 
served explosions, from SGR 1900-1-14 and SGR 1806-20, 
showed similar behavior leading to and following the GFs: 
months before the GF X-ray activity increased, spectrum 
hardened and spindown increased ( |Mereghetti et al. 2005| l. 
In the post-flare period, pulsed fraction and the spin-down 
rate have significantly decreased and the spectrum softened 
(Woods ct al. 2001 Rca ct al. 2005}. AU these effects are 
in agreement with the prediction of twisted magnetosphere 
model l |Thompson et al. 2002} , with the twist increasing be- 
fore the GF and decreasing during the GF brought about 
by reconnection l |Lyutikov 2006| l. Most importantly, obser- 
vations of the December 27 GF from SGR 1806-20 show 
very short rise time ~ 0.25 msec llPalmer et al. 2005t . Such 
short rise times are inconsistent with crustal deformations 
and points to magnetospheric origin of GF jLyutikov 20061 1. 

In analogy with Solar flares magnetic energy to be dis- 
sipated may build-up on long time scales and then be dissi- 
pated on very short times scales, possibly as short as Alfven 
crossing time. For example, a slow plastic motion of the crust 
implants a twist (current) in the magnetosphere on a long 
time scale. At some point a global system of magnetospheric 
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currents and sheared magnetic fields loses equilibrium and 
produces a fiare. This instability may be dynamical {e.g. loss 
of magnetic equilibrium) and/or resistive (this work). 

Though similar in underlying physical process, mag- 
netar and Solar plasmas may difi'er considerably, since 
properties of plasma expected in magnetar magnetospheres 
are very different from the conventional Solar and lab- 
oratory plasmas. The principal difference is that mag- 
netar plasma is strongly (relativistically strongly) mag- 
netized. It is convenient to introduce a parameter Om 
lIKennel fc Coroniti 19841 as the ratio of the magnetic en- 
ergy density ub = /Sn to the total plasma energy density 
(including rest mass!): 



UB 

u„ 



(1) 



Conventionally, in a non-relativistic Solar and laboratory 
plasmas, parameter is small, am ^ 1- On the other 
hand, in magnetar magnetospheres it is expected that am is 
very large. 



~ 10 < (Tm < — 

c \ m„ / il \ rrir. 



10' 



(2) 



(the upper limit here comes from electron-ion plasma den- 
sity equal to the Goldreich- Julian density, the lower limit 
corresponds to plasma density ~ B /{Rnss) at which point 
currents drifting with near the speed of light create toroidal 
magnetic fields of the order of the poloidal). Here lub = 
eB/m^c is cyclotron frequency, B is magnetic field at the 
neutron star surface, Rns is radius of neutron star and Q is 
rotational frequency. 

Large expected value of am (or small 1/am) may be 
used as an expansion parameter in equations of relativistic 
magnetohydrodynamics (MHD). The zero order equations, 
the equations of magnetodynamics (MD), describe the dy- 
namics of magnetic field under the action of magnetic pres- 
sure and tension (IKomissarov 20021 . Since to this order the 
rest mass density and thermodynamic pressure of plasma 
completely vanish, the only effect plasma has on this dynam- 
ics is due to its high conductivity which requires E • B = 
and < B^ . Alternatively, one can write the equations of 
magnetodynamics as Maxwell's equations closed with a suit- 
able Ohm's law that gives the electric current as a function of 
the electromagnetic field only IIGruzinov "19991. This Ohm's 
law is derived from the condition of vanishing Lorentz force 
that explains the alternative name for this system, force- 
free degenerate electrodynamics or simply force-free electro- 
dynamics. 

The limit am — ^ cxd is somewhat reminiscent of the sub- 
sonic incompressible hydrodynamic. In this limit the Alfven 
four- velocity is ua ~ V^am — > oo, so that both cases are 
applicable when (four)-velocities are much smaller than the 
velocity of propagation of disturbances in a medium ( Alfven 
and sound waves correspondingly). In fact, this analogy 
turns out to be much deeper: as we discuss in this paper 
there are deep similarities in governing equations between 
magnetodynamics a small velocity limit of the conventional 
magnetohydrodynamics . 

In order to address energy dissipation and eventually 
particle acceleration, one needs to go beyond the ideal ap- 
proximation. One possible approach is to take into account 
the dissipation of mostly force-free currents due to plasma 



resistivity. Resistivity will result in the decay of currents 
supporting the magnetic field, which would infiuence the 
plasma dynamics. Though plasma resistivity will be of the 
anomalous type, mediated by plasma turbulence and not 
by particle-particle collisions, as a first step we assume that 
plasma resistivity can be represented by macroscopic resis- 
tivity parameter rj. As another simplification we neglect a 
possible back reaction of the heated plasma on the global 
dynamics. This may be justified if the dissipated energy is 
quickly radiated away. 

The principal resistive instability in a conventional, 
non-relativistic plasma is a so-called tearing mode. 
It is one of the principle unstable resistive modes, 
which plays the main role in various TOKAMAK dis- 
charges like sawtooth oscillations and major disruptions 
( Ka domtsev 19751 . unsteady reconnection in Solar fiares 
{ Shivamogg^^^SSIJAschwan deii "200211 and Earth magneto- 
tail ( G alo e^^^^l 9 78 1 . Qualitatively, the most important 
property of the tearing mode is that a current layer of thick- 
ness I may dissipate on time scales much shorter than it the 
resistive time scale ~ /v- In addition, tearing mode may 
be an initial stage of the development of the (steady-state) 
reconnection layers. 

Development of the tearing mode in a relativistic, 
strongly magnetized plasma {am S> 1) is the principal topic 
of this work. Initially, this problem has been considered by 
Lyutikov 120031 . who found that in resistive magnetody- 
namics the current layers are unstable towards formation of 
resistive small-scale current sheets, that's to development of 
tearing mode. He also found that the growth rate of tear- 
ing mode is intermediate between the short Alfven time 
scale Ta ~ l/va (which in a am ^ 1 plasma is the light 
crossing time scale Tc ~ l/c) and a long resistive time scale 
Td'. T ~ {TdTa)^^^. Surprisingly, this is exactly the same ex- 
pression as the one found in the non-relativistic framework 
of incompressible magnetohydrodynamics. In this paper we 
uncover the deep underlying reasons for this coincidence. It 
turns out that slow, resistively-driven evolution of strongly 
magnetized plasma is described by a system of equations 
that is very similar to nonrelativistic MHD. In addition to 
the analytical study we test the theory by means of numer- 
ical simulations. 

In Section|5|we describe the basic equations of resistive 
magnetodynamics. In Section |H] we consider the case of slow 
evolution near equilibrium and show that in this regime the 
MD equations can be reduced to the system that is very 
similar to nonrelativistic MHD. Further reduction is possi- 
ble when the equilibrium is characterized by zero magnetic 
tension. This is shown in Section 2] The analytic theory of 
the tearing instability is reviewed in Section|S] Section|H]de- 
scribes our numerical method and the results of numerical 
simulations are presented in Section [7| 



2 BASIC EQUATIONS 

In any inertial frame of special relativity the dynamics of 
electromagnetic field is described by Maxwell's equations. 
They are the Faraday law 

1 dB 

c dt 

the magnetic Gauss law 



Vx£; = 0, 



(3) 
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V-B = 0, 

the Ampere law 

IdE „ „ 4tt . 

- + VxB = —J, 

c at c 

and the electric Gauss law 
V-E = 4TTg. 



(4) 



(5) 



(6) 



V 



When these equations are supplemented with the Ohm law, 
that relates the electric current with the electric field, the 
system become closed. However, one can also write down 
additional equations that describe the evolution of the en- 
ergy and momentum of the electromagnetic field. These ad- 
ditional equations are very useful for physical interpretation 
of electromagnetic phenomena. 
Energy conservation law: 

de 
di 

where 



+ V-S - 



-Ej, 



(7) 



is the energy density and 
S = —ExB 

is the energy flux density (the Poynting flux) . 
Momentum conservation law; 

^ + V-T = -gE--jxB, 
at c 

where 

is the momentum density and 



(8) 
(9) 

(10) 
(11) 



T = - — (E(g)E + B(g)B) + —{E^+B^)g (12) 

is the Maxwell stress tensor. In the last equation g is the 
metric tensor of Euclidean space. 

In strong magnetic field the conductivity across mag- 
netic field can become highly suppressed (when the Larmor 
frequency becomes much higher than the frequency of par- 
ticle collisions). Moreover, if the magnetic field is so strong 
that one can ignore the inertia of plasma particles then the 
Ohm's law can be written as 

j = (>^^c + ''\\En, (13) 

where iSy is the component of electric field that is paral- 
lel to the magnetic field. Another condition for ea. (ll3ll to 
hold is _E < B. If this condition is not satisfied the cross- 
field conductivity has also to be taken into account, e.g. 
via adding the a±E± term to the right hand side of ea ll3l 
In general, the conductivity o-y, mediated by particle-wave- 
particle interaction, should depend on the Lorentz factor of 
plasma through Lorentz transformation of the effective col- 
lision rates ( |Lyutikov 20 03). As a simplification, applicable 
in cases when no strongly relativistic motion along the mag- 
netic field lines are expected, here we will assume that (T|| is 
a constant macroscopic parameter. 

A characteristic speeds of magnetized plasma is the drift 
speed of charged particles across magnetic field 



ExB 



(14) 



Note that it appears in the Ohm's law 1)13^ where it describes 
the non-conductive contribution to the current density. 

Given the Ohm law 1131 one can write the energy and 
momentum conservation laws as 



+ V-S 



(15) 



(16) 



VxE^O, 



(17) 



de 
di 

f + V.T = -,B„ 

In the limit of infinite conductivity, (7|| oo, the electric 
current — o-||_E|| must remain finite and this ensures oc 
0. In this limit the dynamics of the electromagnetic 
field can be described by the following closed set of equations 
IKomissarov 20021 

1 dB 
c~dt 

V-B = 0, (18) 
|^ + V-S = 0, (19) 

^ + V-T = 0. (20) 

In (Komissarov 2002| this system of equations was derived 
from the system of ideal relativistic MHD in the limit of van- 
ishing rest mass density and pressure of matter. This way 
of derivation suggests to call this system magnetodynam- 
ics (MD), the name that is obtained from magnetohydrody- 
namics (MHD) via excluding its hydro part. Alternatively, 
it could be called force-free degenerate electrodynamics. This 
name originates from the early attempts to construct steady- 
state models of magnetospheres satisfying the degeneracy 
condition, 

E B = 0, 
and the force-free condition. 



gE+-jxB 

c 



0, 



e.g. IThorne fc Macdonald 198211 . However, it is longer and 
makes an unwarranted emphasis on the electric component 
of the electromagnetic field. Indeed, similarly to ideal MHD 
where the electric field vanishes in the fiuid frame, in ideal 
MD the electric field vanishes in the frame moving with the 
drift velocity. Thus, the electric component of the field can 
still be considered as a secondary one that is induced by 
the motion of magnetized plasma. In fact the electric field 
vector can be eliminated from the set of dependent variables 
of equations II17I20I I via 



E 



1 



VxB, 



1 



T = 



:{VxB) ® (VxB) 



+ 



(21) 
(22) 

(23) 
(24) 
(25) 
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1 + 



Moreover, in the dissipative regime the electromagnetic field 
can no longer be described as force- free either (see ea ll6l '). 



3 THE QUASI-EQUILIBRIUM 
APPROXIMATION. 

Here we will derive equations that describe slow evolution 
of strong electromagnetic field near the state of force-free 
equilibrium. There are two basic characteristic time scales 
in such problems, the light crossing time scale 



Tc = I/C 

and the diffusion time scale 



Td 



(26) 



(27) 



where rj = /4na\\ is the resistivity and I is the character- 
istic length scale of the problem. Since in the limit of mag- 
netodynamics the hyperbolic waves (fast and Alfven waves) 
propagate with the speed of light then Tc is the typical time 
of establishing dynamical equilibrium. Td gives the time scale 
for the diffusion of magnetic field due to finite resistivity. For 
systems close to a static equilibrium the characteristic time 
scale of their global evolution, which we will denote as r, is 
much larger than Tc. This gives us the small parameter 



/ic — Tc/t <C 1. 



(28) 



Diffusion is usually only important on small scales where 
there may exist large gradients of physical quantities. There- 
fore we will also assume that r is much smaller than Td and 
this gives us another small parameter 



fj-d ^T/Td 1. 



(29) 



This automatically ensures that the relativistic Lundquist 
number 



Lu = V2Td/Tc = V2lc/-n > 1 



(30) 



(\/2 is introduced in order to simplify the expressions of 
Sec|S]) Moreover, from the Faraday equation one immedi- 
ately obtains 



(31) 



where £ and B are the characteristic scales of the electric 
and the magnetic fields respectively. Thus, the electric field 
is much weaker than the magnetic field and this results in 
the drift speed much smaller than the speed of light (see 

The Ohm law 11311 leads to the following dimensionless 
form of the Ampere equation 



2dE „ „ 2 ExB 



where 



5i = (47r)"L„-^. 



Thus, the Ampere law may be reduced to 



(32) 



(33) 



(34) 



and the characteristic scale for the parallel component of the 
electric field is 

= (47r)-^L-'S < B. (35) 

Using this equation and the expression (1141 for the drift 
velocity we can write 

E = --VxB + -^VxB. (36) 
c 47r(T|| 

Substitution of this result into the Faraday equation leads to 
the exactly the same advection-diffusion equation for mag- 
netic field as in the nonrelativistic MHD 
dB 



dt 



Vx {VxB)-rj\/^B = 0. 



(37) 



The dimensionless form of the energy equation H15^ 
reads 

^ + V-S^~S2El (38) 
where 

- « 1- (39) 

Thus, the dissipative term may be omitted and the energy 
equation reduces to its ideal form 

de 



dt 



(40) 



The dimensionless form of the momentum equation I16II is 

(41) 



(42) 



^il^ + v■T = ^53pE^^, 

where 

% = 7^ « 1- 

Thus, the resistive term in the momentum equation can also 
be omitted and we write 

f + V.T = 0. 

ot 

Given condition 13H we may write 



8n 



and 



1 B 
B^B + g — . 

47r 87r 



(43) 



(44) 



(45) 



Then we combine eq.||lj, eq.J^ZJ, iQOli a^nd H43^ into a closed 
system of equations very similar to nonrelativistic resistive 
MHD 



V B = 0, 

^ - V X (V X B) - rjV^B = 0. 
| + V(2pV).0. 

In these equations 



87rc2 



(46) 
(47) 

(48) 
(49) 

(50) 
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is the mass density of the electromagnetic field. There are 
however several differences. First, there is the factor of 2 in 
the second term of the "continuity equation" (1481 . Second, 
there is no V(pV ® V) term in the momentum equation. 
Finally, there is no "energy equation". To some extent this 
similarity has been noticed in ()Gmzinov 1999^ . 



4 THE INCOMPRESSIBLE LIMIT. 

This system can be reduced even further if the equilibrium is 
supported predominantly by the magnetic pressure. Indeed, 
the dimensionless equations of momentum and continuity 
include small parameter e = jj.^ 



dp 
at 



+ 2V(pV) =0, 



(52) 



Using the method of perturbation we expand the unknowns 
in powers of e 

B = Bo + eBi, y = Fo + eFi, p = po + epi, 

where 

(Bo-V)Bo = 0. (53) 
To the zero order we have 



87rc2 



const 



(54) 
(55) 



and 

V-Vo = 0, 

whereas the first order gives us 

Po^-^(Bi-VBo+Bo-VBi) + V:?i-^=0. (56) 
This result shows us that eg. 15111 difi^ers from 

only by terms of order and we may approximate equation 
lUnt by 

that implies the vorticity equation 

dV xV 1 

^^^-4;^^^(^-^)^ = °- 

Thus, we have arrived to the following system of equations 
V-B = 0, (60) 
V-F = 0, (61) 



^ - V X (y X B) - TjV^B = 0, 
P^-ivx(B.V)B^O. 



dt 



4-K 



(62) 



(63) 



This system of equations is very similar to the one of incom- 
pressible MHD. The only difference is that equation I58II 
involves the Eulerian time derivative whereas the vorticity 
equation of incompressible MHD involves the Lagrangian 



time derivative (this, however, has no effect on the pertur- 
bation equations of the linear stability theory). One may 
argue that the limit of incompressible magnetodynamics is 
applicable only to a very limited range of problems as in 
general case the magnetic tension cannot be ignored and as 
the result one cannot assume that p is constant. We only 
point out that exactly the same argument applies to incom- 
pressible MHD unless the plasma magnetization is very low 

l |Bisl<:amp 1993t . 

Both in the nonrelativistic case llPriest fc Forbes 2000t 
and in the case of magnetically dominated relativistic 
plasma l |Lyutikov 2003^ the most basic equilibrium cur- 
rent sheet configuration involves 1-dimensional current sheet 
with slab symmetry. Within this current sheet the mag- 
netic field gradually rotates so that its pressure remains 
constant and its tension vanishes. Studying the stability 
of the relativistic current sheet Lvutikov ll2UU3t derived the 
same growth rates for the tearing mode as in the nonrel- 
ativistic theory. Our analysis of the equations of resistive 
magnetodynamics explains this result as the direct conse- 
quence of the similarity between magnetodynamics and non- 
relativistic MHD. In the following section we study the tear- 
ing mode in some details. 



5 TEARING MODE INSTABILITY 

Consider the stationary versions of eas. 1621631 

Vx (VxB)-r?V^B = 0, (64) 



—V X (B-V)B = 0. 

47r 



(65) 



Low 119731 found the following one-dimensional solution for 
these equations 



B = Bota,n]i{y/l)ix i Bosech{y/l)iz, 
V = ~{v/l)tanh{y/l)iy. 



(66) 
(67) 



It describes a current sheet in which the magnetic field ro- 
tates exactly by n without changing its magnitude. The ro- 
tation continues from y — — oo to y = +oo but most of 
it occurs within y G [—1,1]- This allows us to describe I 
as the thickness of current sheet. The total electric field , 
E = — V xB + rj'V X B, corresponding to this solution is 



(68) 



It turns out that the electric and magnetic fields given 
by equations 1)67168^ also satisfy the full set of stationary 
Maxwell's equations supplemented with Ohm's law 1131 . (As 
a matter of fact this provides us with a very good test prob- 
lem for our numerical code.) 

Consider a perturbation of this equilibrium state of the 

form 

b(x, y, t) = b{y) exp{i{ka:X + kyy) + wt), 

v{x,y,t) = v{y) exp{i{kxX + kyy) + wt). 

Equation 1621 and curl of equation 1631 give us the following 
linearized equations for the the y-components of the pertur- 
bation 
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why = ivy{k-B) +7]% - k by), 



w{iy 



i(k-B) 



-b, 



k B 
k B 



+ {by 



(69) 



(70) 



These equations are the same as equations 6.6 and 6.7 in 
Priest & Forbes |J2000| where the tearing mode instability 
is studied in the non-relativistic framework of incompress- 
ible MHD (Priest & Forbes I2U001 reproduced the original 
analysis by Furth et al. 11963(1 in a somewhat more concise 
and transparent way.) Thanks to this, the nonrelativistic 
theory and the relativistic theory should deliver exactly the 
same results. At this point we could simply refer the reader 
to the studies cited above but for the sake of completeness 
we briefly outline the analysis of Priest & Forbes II2UUUII . 
clarifying few important issues along the way. Moreover, in 
our version we will use the steady state solution given by 
eas. 16716711 instead of the step function 

r +Bo for x>l, 
B:r = < Box for !a;| < Z, 
[ —-Bo for X < —I. 

used by the cited authors for the sake of simplicity. 

Equation 17U1 is singular at the point where the wave 
vector is normal to the field direction in the equilibrium 
solution, that is where k - B = 0. At this location there 
develops a thin resistive sub-layer where the second order 
derivatives become large just like in the classical boundary 
layer problem (e.g. Bush, 1992). With some loss of generality 
we will assume in what follows that fcz = 0, so that k ■ 
B = kBx- This puts the sub-layer right in the middle of 
the current sheet and makes the problem symmetric with 
respect to the plane y = 0. 

Using / as a unit of length, Td — /rj as a. unit of time, 
and Bo as a unit of magnetic field strength we can write the 
dimensionless version of equations 1691761 as 



Wby = -Vyf + {by -k by), 



5{vy -k Vy) ^ f 
where 



-by^ + (by-k%) 



(71) 



(72) 



b = —— , k = kl, w — 



kl^ 



Bo' 



, Vy = -tVy , 



and 



/(y) = ^ = tanhy, y = y/l, 
Bo 



Luk"^ 



is a small parameter. Expanding by and Vy in powers of 5 
we find the the leading terms satisfy 



by — k by + 2(1 — tanh y)by = 0, 
^(2tanh^y-2- w). 



(73) 

The solution to equation 17311 that satisfies the boundary 
condition by{+oo) = is relatively simple 

fc/2 



6o(tanhy + k) / I — tanhj/ 
k V ^ + tanh y 



(75) 



Given by one can find Vy from eg. 17111 . Once by and Vy are 
known one can find bx and Vx from the conditions 

V-b = 0, V-w = 0. (76) 

From eg. 17511 we find that 6(0) — bo and then eg. l7H 
shows that Vy diverges at y = 0. This indicates that around 
y = there exists a sub-layer where the second order deriva- 
tives of Vy in equation 16911 become large and cannot be 
dropped. Thus, solution H75|l applies only to the outside of 
this sub- layer and is called the outer solution. In fact, it holds 
only for y > 0. The solution that applies for both sides of 
the sub- layer is 



bo (tanh \ y\ + k) 



1 — tanh \y\ 
1 -I- tanh \y\ 



fc/2 



(77) 



This solution is continuous ai y — but its first derivative 
changes sign at this point so that 
n 0+ 



b' 



= 2- 



1-fc^ 



k 



(78) 



These derivatives are to be used for matching the outer so- 
lution with the inner solution that describes the interior of 
the resistive sub-layer. 

To find the structure of the resistive sub-layer, or the 
inner solution, we need to introduce new stretched variable 
^ = y/S^ , where the power of S is to be found using the prin- 
ciple of least degeneracy \Van Dyke 1975| |. Once p is known 
we will find thickness of the sub-layer as e = S^. Here we 
may assume that inside the sublayer y <C 1 and use the ap- 
proximation f — y. Combining equations 17111 and 1721 we 
find 



w y 



k^v. 



y , y - 

W 



(79) 



Substituting this result into eg. 1711 we end up with a 4th 
order ODE for Vy. When using the stretched variable this 
equation reads 

^2^(4) _ 2^(3)^ + 25^^-^v'^^'> 

-S^'^-^i'^ik^ + 2w)v = A{i, v) 



(80) 



where A(^,v) includes all the term with positive power of 
5. One can see that when Vy is expanded in powers of 5 
the leading order equation will retain the largest number of 
terms if p = 1/4. Thus the dimensionless thickness of the 
sub- layer is 



, 1/4 



f;(2) 



f— V 

\Llk-^) 

When using the stretch variable eg. 1711 reads 

\W + wYby - iVy =0. 



(81) 



(82) 



To the second order in e we can ignore the last term in this 
equation and write 



f,(2) 



\k^ + W)by = 0. 



(83) 



The solution of this equation that is symmetric with respect 
to the ^ — 5 transformation is 



by — bo cosh(-\/ fc2 -f wt£). 

this gives us the jump in derivatives of by across the sub- 
layer 
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A' 



6y 




Uy 




by 


e 

y=-€ 


. by . 





2(fc^ + w)e. 



(84) 



Here by is the derivative with respect to y. 

Matching A's derived from the outer and the inner so- 
lutions, eg. 18411 and eg. 1781 respectively, allows us to find 
the growth rate for the mode with wavenumber k: 



(k -\- w 



1/4 



(85) 



One can see that only modes with wave numbers 

k < 1 (86) 

grow in amplitude. For fc ^ 1 dispersion relation simplifies 
resulting in the well know equation 
/ r X 2/5 

According to the last result the growth rate increases 
indefinitely with the wavelength. However, the e-folding time 
of the growing mode 

f 



cannot be shorter than the resistive time scale of the sub- 
layer 



T^V 2 

This condition reads us 



(89) 



(90) 



For smaller k the growth is reduced and thus the maximum 
growth rate should be of order 

- rl/2 
w — wz, = L,; 



(91) 
(92) 



that corresponds to the dimensional time-scale 

• f ^l/2 

T = (TcTd) 

6 NUMERICAL METHOD 

To verify the predictions of the theory of the tearing mode 
instability in relativistic magnetically dominated plasma we 
have carried out a series of numerical simulations. The equa- 
tions that are solved numerically are the Maxwell equations 
13I6II with the Ohm law 

E X 13 

3 = Q^:^c + cj\\E\\+o^E^. (93) 

This differs from eq. 11311 by the presence of the conductivity 
current across the magnetic field. However, we set a± to 
when < and thus eg. 19311 reduces to eg. 11311 during 
the linear and the beginning of the non-linear phases of the 
tearing instability. However, at some point of the nonlin- 
ear phase a region develops where the electric field becomes 
stronger than the magnetic one thus allowing significant 
cross-field conductivity. We expect the particle pressure and 
inertia to become important in such regions resulting in a 
different form of the Ohm law. In fact, the MHD approxima- 
tion may provide a more suitable mathematical framework 
there. However, in this study we are more focused on the lin- 
ear theory and for this reason the prescription 1931 suffice. 



The numerical scheme used for simulations is a derivative of 
the general relativistic scheme for resistive magnetodynam- 
ics described in Komissarov 120U4t . Here, we give its brief 
description. 

From eq.© one finds 



at(V-B) = 0. 



(94) 



This well known result shows that it is sufficient to enforce 
the divergence free condition Q only for the initial solu- 
tion and it will be satisfied automatically at any other t. 
Unfortunately, straightforward application of many numer- 
ical schemes perfectly suitable for other hyperbolic systems 
of conservation laws fails to deliver a good result for elec- 
trodynamics and MHD simply because their discrete equa- 
tions are not consistent with any discrete analogue of 19411 . 
In particular, this applies to the method of Godunov which 
has many beneficial properties and is currently considered 
as generally superior to many other numerical schemes for 
hyperbolic systems. There have been many attempts to find 
a cure for this "div-B problem" ( see the review in Dedner 
et al.,2002.) 

One of the ways to handle this problem involves con- 
struction of a somewhat different system of differential equa- 
tions, the "augmented system", where the divergence free 
condition @ is no longer included and V-B may be trans- 
ported and/or dissipated like other dynamical variables. The 
idea is not to enforce the divergence free condition exactly 
but to promote a natural evolution of the system toward 
a divergence free state. Provided the augmented system is 
hyperbolic it can be solved numerically using the Godunov 
method (Munz 2000 Dedner ct al. 2002 ). 

Following this idea we modify eq.Q as follows 

+ Vx£; + V* = (95) 

c at 

where 'I' is the scalar field which is called the pseudopotential. 
It's evolution is given by 

1 d<f 



dt 



+ V-B + Kc* = 0. 



(96) 



From these one finds that both ^' and V • B satisfy the 
telegraph equation 

1 a^* 



9£ 
'dt 



0. 



+ K — hV^(V-B) =0. 



(97) 



(98) 



c2 dt^ dt 

The resultant system includes two vector equations 
195I5I I. and one scalar conservation law, equation 1961 . All 
these evolution equations can be written as conservation 
laws. In the form of an abstract vector equation they read 



dt dx^ 
where 

Q = {'i',B\E') 

are the conserved variables 

= c {B\e''''Eu + *<7'^ -e'^'=Bfc) 
are the corresponding hyperbolic fluxes, 
S = («:c^*,0\47rj'). 



(99) 
(100) 
(101) 
(102) 
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In this study we employ Cartesian coordinates so that 1) the 
metric tensor of space g*-' equals to Kronecker's delta (5'-' , 2) 
its determinant 7=1, and 3) the Levi-Civita alternating 
tensor e*-** reduces to the Levi-Civita symbol with e^^^ — 1. 

Numerical integration of 19911 is carried out using a sec- 
ond order Codunov method based in combination with the 
time-step splitting technique by Strang l| 19(i8|l for handling 
the source term. Namely, each time-step — > -I- /S.t is 
split into three sub-steps. During the first sub-step the nu- 
merical solution is advanced in time by At/2 via integrating 
the truncated system 



In fact we split the source term into two parts 
Sa = (kc^*, O'ATifc) and Sb = (0, Q\ Anu), 



(103) 



(104) 



where = (^\\E\\ + a±E± is the conductivity current and 
= QcExB/B is the drift current, and split this substep 
as well. First we account for Sa only. This source term is po- 
tentially stifle but its simple form allows exact analytical inte- 
gration of ea. 110311 thus reducing the stability constraints on 
the timestep. Then we integrate eg. 110311 with source term 
Sb numerically using the method of Newton. During the sec- 
ond sub-step the resultant solution is advanced in time by 
At via integrating the truncated system 



dt 



= 



(105) 



using the second order Codunov scheme. The third sub-step 
is a repeat of the first one but now it is the solution found 
by the end of the second sub-step that is used as the initial 
solution for eg. 110311 . 

In our Codunov scheme the numerical fiuxes through 
the cell interfaces are found using the exact solution to the 
interface Riemann problems. The construction of this exact 
solver is simplified by the fact that system 11051 is linear. 
Its ID-version can be written as 



dt dx 







(106) 



with the Jacobean matrix 



A = c 



The eigenvalues of this matrix 



- 


1 

















1 
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-1 















1^(1) = M(2) = M(3) = c; 



M(4) = 0; 



M(5) — M(6) — M(7) — 



provide the wavespeeds of hyperbolic waves. Other proper- 
ties of these waves are given by the eigenvectors of A. The 
right eigenvectors are 



r(i) = (1,1,0,0,0,0,0) 

r(2) = (0,0,-1,0,0,0,1) 

r(3) = (0,0,0,1,0,1,0) 

r(4) = (0,0,0,0,1,0,0) 

r(5) = (-1,1,0,0,0,0,0) 

r(6) = (0,0,1,0,0,0,1) 

r(7) = (0,0,0,-1,0,1,0) 



(107) 



Since A is symmetric its left eigenvectors, coincide with 
the corresponding right eigenvectors: 



hi 



(108) 



It is easy to see that solutions 2,3,6, and 7 are the usual 
electromagnetic waves. Solution 4 simply refiects the fact 
that all waves of vacuum electrodynamics are transverse and 
any discontinuity in the normal component of electric field is 
due to a surface electric charge distribution. Solutions 1 and 
5 describe new waves that do not exist in electrodynamics; 
they transport ^ and V-S. 

The solution to the Riemann problem with the left and 
the right states, and Q(r), and the interface speed (3^ is 



Q = 2(0 + 



i=l,3 



(109) 



where the wave strengths 

(2m 



2(o)-r(o 



(i) ■ ' (i) 



Notice, that the 4th wave is ignored. The x-component of 
electric field is set to be the mean value of the left and the 
right states 



E"" = 0.5{E(i) + E\ 



(r)) 



(110) 



The entire scheme is second order accurate in space and 
time provided /S.t ^r\lc. The only stability constraint on the 
time step comes from the Currant condition At < /S.x/c for 
eg. 11051 . However, the accuracy considerations may require 
a smaller time step. For example, our ID test simulations 
of the equilibrium current sheet 167I67II have shown that if 
At = rj/c then the relative error in E^ is around 10%. 



7 COMPUTER SIMULATIONS 

These are 2D simulations with slab symmetry in the z- 
direction (that is df/dz — for any function /). The utilized 
units are such that / = 0.1, c = 1, and Bo = 1. Following the 
setup of Sec|3the initial equilibrium current sheet is normal 
to the y-axis and its symmetry plane is j/ = 0. The compu- 
tational domain is [— xc+xq] x [— yo,+?yo] with xo = A/2, 
the half of the perturbation wavelength, and yo = 2.0. On 
the boundaries x = ±10 we impose the periodic boundary 
conditions and on the boundaries y = ±yo we impose the 
zero-gradient conditions. In order to get accurate solutions 
we have to have the resistive sub-layer well resolved. On the 
other hand, far away from the sub-layer the solution has no 
fine structure and does not require high resolution. This sug- 
gests to use a grid with variable resolution in the y-direction. 
We adopted the exponential rule 



aAy^ 
a-^Ay, 



for 
for 



<0, 
>0, 
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Figure 1. The evolution of the current sheet with k ■ 



: 0.314 {kLi^* 



1.08). Bz (color image) and the magnetic field lines at t = (top 



left), t = 5 (top right), t = 10 (bottom left), and t = 16 (bottom right). 



with the resolution in the symmetry plane being 10 times 
higher than that at the domain boundary, Aj/(0) = 
0.1Ay{yo). The resolution in the x-direction was uniform 
with Ax = VoJAy{yo). 

We have tried two different models for the initial solu- 
tion. In the first model the perturbation of the equilibrium 
solution H67I68|I has the form 



(0,6osin(7ra;/2;o),0), 



(0,0,0) 



with bo = 10""^. Since this is not the normal mode of the 
tearing instability the corresponding numerical solutions ex- 
hibited initial settling period of order 5t ~ 1 before reaching 
the stage of exponential growth. 

In the second model we tried to set up the normal mode. 



In fact, we used the the outer solution 1751741 with bo = 10~^ 
in order to introduce the perturbation. Since, the electric 
field given by this solution diverges at yo we resorted to lin- 
ear interpolation in order to continue it within the resistive 
sub- layer. Perhaps this explains why the solution still exhib- 
ited approximately the same settling time as in first model. 
Moreover, we have not found any significant differences be- 
tween the solutions corresponding to the both models except 
from their behavior during the settling. 

We have considered only two values for the resistivity, 
77 = 10~^ 77 = 10"*, corresponding to the Lundquist number 
Lu = 1.4 X 10^ and L„ = 1.4 x 10^ respectively. Figure Q 
shows the evolution of the current sheet for ri = 10~^ and 
k = 0.314. As the perturbation grows the current sheet grad- 
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Figure 2. Multiple magnetic islands at the nonlinear stage of the 
model with k = 0.0393 (kL]^"^ = 0.135, r) = IQ'^). 

ually thins out in the middle of the domain and thickens at 
the x-boundaries. Eventually this results in the development 
of two large magnetic islands. By the end of the simulations a 
third much smaller island forms right in the middle. In fact, 
such secondary islands are typical and more pronounced for 
longer wavelengths. Figure |21 shows the final solution for the 
model with k = 0.039. In this case there are free relatively 
large secondary islands and few smaller secondary islands in 
between them. These secondary islands originate in the fol- 
lowing order. The larger secondaries develop when the two 
primaries suck in a large fraction of the current lines so that 
the residual current sheet becomes rather thin. These first 
secondary islands also suck in the current lines and this leads 
to further reduction of the thickness of the residual. Then 
new smaller islands appear between secondary island of the 
first generation and so on. This behavior seems to be caused 
by the reduction of the growth time for the tearing mode in 
thin current sheets. 

As the thickness of the residual current sheet reduces, 
V X _B goes up and so do the electric field and electric cur- 
rent density in the sheet. At some point the electric field, 
which is predominantly perpendicular to the magnetic field, 
reaches the same strength as the magnetic one and the drift 
speed speed reaches the speed of light. Eventually the cur- 
rent sheet collapses to the size of a computational cell and its 
electric field becomes even stronger than the magnetic field. 
This shows that the structure of the current sheet can no 
longer be described within our model and other factors like 
the thermodynamic pressure of plasma heated to high tem- 
peratures in the current sheet has to be taken into account. 
Prior to this point the quasi-equilibrium of the current sheet 
is basically supported by the magnetic pressure alone, the 
pressure of predominately x-directed magnetic field outside 
of the current sheet being balanced by the pressure of pre- 
dominately z-directed magnetic field inside of it. 

Figure |3] shows how the maximum value of By grows 
with time for different wavelengths of the perturbation in 
the case of 77 = 10"^ One can clearly see that after the 
short initial settling period the growth of By becomes expo- 
nential. At this phase the primary islands dominate in the 
solution. This phase terminates when the nonlinear effects 




5 10 15 20 25 30 

t 



Figure 3. Growth of the perturbation for various wavenumbers 
{rj = 10^^). BJ'"" is the maximum value of By on the grid. 

become significant or when the smaller secondary islands 
overtake the primary ones (see the curves for k = 0.0393, 
and 0.0196). There are two reasons for the faster growth of 
the secondary islands. Firstly, at the point of their appear- 
ance the thickness of the current sheet has been significantly 
reduced. Secondly, both the primary and the secondary is- 
lands correspond to wavenumbers fc < fc*. In this regime 
(see figure the perturbations of shorter wavelengths grow 
faster. 

Figure 3] shows the dispersion relations w = w{k) that 
accounts for the growth rates of the primary islands. One 
can see that the dispersion curves have maximum at around 
k* — Lu^^'^ that agrees very well with the analytical results 

presented in Sec|5] The maximum growth rate is also close 

— 1/2 

to the predicted wg. = LJ . Moreover, the results clearly 
indicated that modes with k > kc — 1 are indeed stable (see 
ea l85ll . The reason why the results for t] = 10~* give the 
maximum growth rate somewhat closer to the theoretical 
value than those for ij = may have to do with the fact 
that the separation between the cut-off wavenumber, kc, and 
k* becomes too short for relatively low Lundquist numbers. 



8 APPLICATION TO MAGNETAR FLARES 

One of the main unknowns in our calculations is the value 
of the resistivity rj, which should be calculated from the par- 
ticle kinetics, but instead was introduced as a macroscopic 
property of plasma - a common approach in continuous me- 
chanics. Due to very short radiative decay times in magnetar 
magnetospheres the particles are bound to move only along 
the field lines, this limits considerably a number of possible 
resonant wave-particle interactions that can lead to develop- 
ment of plasma turbulence. The main remaining options are 
the Langmuir turbulence, which in relativistic plasma devel- 
ops on a typical scale of electron skip depth, 5e ~ c/up^e, 
and ion sound turbulence, which in relativistic plasma de- 
velops on an ion skip depth Si ~ c/ujp^i (t<Jp,e and ujp^i are the 
electron and ion plasma plasma frequencies). They are dif- 
ferent by a square root of the ratio of electron to ion masses 
5e/5i ~ (me/jTii)^''^. A fuUy developed turbulence with a 
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Figure 4. Growth rate as a function of the wavelength. Dashed 
line: rj = 10~^ solid line: rj = 10~*. 



typical velocity c and typical scale 5 would produce a resis- 
tivity Tj ~ c5. Thus we can write the growth time of tearing 
mode as 



Notice, that the result does not depend very much on which 
of the skin depths we use. The difference is in the fac- 
tor {m^lm,i)^l^ ~ few. In current-carrying magnetospheres 
of magnetars the plasma frequency can be estimated as 
ojp = X^Jujbt /c (A is a parametrization constant and r is 
the radial distance from the neutron star center) giving a 
growth time of the tearing mode 

r = 2{B/BQf'\l/RNs)''*^ msec, (112) 

where Rns is the neutron star radius and Bq is the quantum 
magnetic field. We also note that the expected Lundquist 
number is L„ ~ ~ 5 x IG^iB/BQf^il/RNsf'^X'^'^, 
which is close to the values used in our numerical simula- 
tions. 

Although the connection between the growth time of 
the tearing mode and the observed flare rise time may not be 
that straightforward, the fact that the estimate 111211 is close 
to a typical rise time of magnetar flares ( IGogiis et al. 1999| l 
is encouraging. Keeping in mind the uncertainties in the 
estimates of r] and I we conclude that our model is consistent 
with observations. 



9 CONCLUSIONS 

The original idea of this study was to verify the surpris- 
ing analytical results of jLyutikov 2003| l on the tearing in- 
stability in magnetically dominated relativistic plasma by 
means of direct numerical simulations. We have carried out 
such simulations and they fully confirmed Lyutikov's conclu- 
sion on the growth rate in magnetically dominated regime 
- it indeed coincides with that found earlier for the tear- 
ing instability in non-relativistic incompressible MHD, the 
shortest growth time being equal to the geometric mean 
of the Alfven and resistive timescalcs (Furth et al. 19631 



IPriest fc Forbes 200011 . Moreover, like in the non-relativistic 
case we observed formation of magnetic islands in the non- 
linear phase of the instability. This is accompanied by the 
appearance of noticeable O-type and X-type neutral points 
(or lines). The development of an X-point may be consid- 
ered as a first step in setting up of reconnection layer where 
dissipated magnetic energy and magnetic flux are constantly 
replenished due to plasma inflow ( Syrovatskii 1 98T| . 

Our results support the possibility of magnetar 
flares being magnetically driven reconnection-like events 
that occur in magnetar magnetospheres l |Lyutikov 2003| 
[Lyutikov 2006^ in a similar fashion to flares and coronal 
mass ejections of Sun's magnetosphere. More complicated 
magnetic field geometries, e.g. those with nonvanishing mag- 
netic tension, should also be subject to tearing instability on 
time scales intermediate between the resistive and Alfven 
time scales llPriest fc Forbes 200011 . The growth rate typi- 
cally scales as r ~ ''""t"^"", where < a < 1 is some coeffi- 
cient (q = 1/2 in our case). In the case of magnetar flares the 
observed flare rise-times range from fractions of a millisec- 
ond to ~ ten milliseconds, which is indeed intermediate be- 
tween the Alfven crossing time, ~ 0.03 msec, and the resis- 
tive time, ~ seconds, of the magnetosphere l |Lyutikov 20031 1. 

The short timescale of tearing mode is achieved entirely 
through the formation of a very narrow resistive sub-layer 
with very short diffusion time scale. However, it is also ex- 
pected that the resistivity itself may be enhanced inside 
such sub-layers due to the development of plasma turbu- 
lence (anomalous resistivity). Elucidating the properties of 
anomalous resistivity in such plasma is an important next 
step. Moreover, the collapse of current sheets and the even- 
tual development in our simulations of the drift speeds ex- 
ceeding the speed of light clearly indicates a breakdown of 
the MD approximation. In order to move forward one needs 
to put back the dynamical effects of matter which are bound 
to become important in this singular domain. The simplest 
way of doing this is to return to the model of relativistic 
MHD like in Lyutikov fc Uzdensky 1200311 and Lyubarsky 
(2005). 

The coincidence of the growth rates for the tearing in- 
stability in incompressible MHD and MD clearly suggests 
some hidden similarity between the evolution equations of 
these two systems. We have carried out further analysis of 
resistive MD and discovered that in the quasi-equilibrium 
limit, that is characterized by low drift speeds, its equations 
do indeed become similar to those of non-relativistic MHD. 
In particular, the mass-energy density of the magnetic field 
takes on the role the mass density of classic MHD whereas 
the drift velocity plays the role of the fluid velocity. So the 
dynamics of the system can be described as a flow of mag- 
netic mass-energy under the action of magnetic pressure and 
tension. Moreover, the fact that the magnetic tension of the 
equilibrium current sheet subjected to tearing instability is 
zero allows further reduction of the MD equations which 
leads to almost the same system as incompressible MHD. 
When, this system is utilized in the analysis of the tearing 
instability it generates exactly the same equations for the 
perturbation as in incompressible MHD. Thus, our analysis 
provides perfect mathematical explanation for the surprising 
results in Lyutikov (2003J S'ld improves our understanding 
of the dynamics of magnetically-dominated plasma. 

Finally, we seem to have come up with a more suit- 
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able name for the system of equations describing the dy- 
namics of relativistic magnetically dominated plasma. In- 
stead of the currently used force-free degenerate electrody- 
namics (FFDE) or simply force-free electrodynamics (FFE) 
we propose to use magnetodynamics (MD). Let us sum- 
marize our arguments in favour of this new name. The 
current name originates from the early studies of station- 
ary magnetospheres of neutron stars and black holes. Be- 
sides the steady-state Maxwell's equations, the key equa- 
tions used in those studies were the degeneracy condition, 
E ■ B — 0, and the condition of vanishing Lorentz force. 
The full system of the corresponding time-dependent equa- 
tions was not known at the time and so were not known its 
connection with the relativistic MHD and the properties of 
it's waves. All these have been discovered only very recently 
llUchida 19971 in ruzmov "19991 IRomissarov 2002ll . In Komis- 
sarov II2UU2II the dynamical equations of relativistic magnet- 
ically dominated plasma, equations II17I2UII . have been de- 
rived from the system of ideal relativistic MHD in the limit 
of vanishing rest mass density and pressure of matter. This 
way of derivation immediately suggests simply to remove the 
/lyrfro-component from the word magneto-hydro-dynamics, 
which obviously results in the name we propose here. Like 
in ideal MHD, in ideal MD the electric field also vanishes in 
the the fluid frame, or to be more precise in the frame mov- 
ing with the drift velocity. An inertial observer moving with 
this velocity would experience a pure magnetic field. Thus, it 
is desirable to retain the emphasis on the magnetic compo- 
nent of the electromagnetic field, that is present in the name 
of MHD. The key wave of MHD, the Alfven wave, is also 
present in MD. It propagates along the magnetic field lines 
with the speed of light ( Komissa rov 2002t . Like in MHD the 
electric field vector can be eliminated from the set of depen- 
dent variables of equations 117I20II (see Sec|5J. Finally, as 
we have shown in this paper one can qualitatively describe 
the evolution of magnetically dominated plasma as a flow of 
magnetic energy under the action of Maxwell's stress. Ap- 
parently, we are not the only ones who are not particularly 
happy with the name force-free electrodynamics. Recently, 
A.Spitkovsky proposed to use the name force-free MHD in- 
stead jSpitkovsky 2006 ) . This name is somewhat better but 
it is still not completely satisfactory. First of all, matter and 
the electromagnetic field appear in this name on equal terms, 
contrary to the relevant physical conditions. Secondly, the 
name MHD implies a somewhat different set of evolution 
equations. Thirdly, as we have seen, the resistive case can 
no longer be described as force-free. Finally, the name mag- 
netodynamics is made of just one word that makes it more 
esthetically pleasing. 
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